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Introduction 


This  report  covers  the  two  major  activities  investigated  this  year. 

(1)  Preprocessing,  feature  extraction  and  classification  algorithms 
have  been  applied  to  classify  data  from  the  NRL  data-base.  In¬ 
variant  moments  and  normalized  Fourier  descriptors  were  the  features 
used  in  these  classification  experiments,  which  show  that  the 
Fourier  descriptors  are  better  suited  to  classification  of  air¬ 
craft  in  turbulent  atmosphere.  Weak  points  of  the  algorithms  are 
discussed  and  directions  for  further  research  are  indicated. 

(2)  Post-processing  of  image  data  to  remove  the  effects  of 
atmospheric  turbulence  was  also  studied.  Atmospheric  turbulence 
was  simulated,  by  generation  of  random  phase  components,  and 
applied  to  an  idealized  model  of  the  drone  in  the  NRL  data-base. 

An  algorithm  for  atmospheric  deblurring,  known  as  shi  f  t-and-add , 
was  used.  The  shi ft-and-add  algorithm  produced  images  free  of 
turbulence  affects  but  degraded  by  blur.  Image  restoration  was 
applied  to  the  processed  images  and  definite  improvements  in 
image  quality  were  obtained. 
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Data  Base  Analysis 

NRl  has  supplied  a  data  base  of  FLIR  imagery  of  a  station¬ 
ary  drone  acquired  at  White  Sands  Missile  Range  (see  Final 
Technical  Report  No.  DIAL-82-010,  October  30,  1982).  Twenty 
images  from  each  of  the  eleven  distinct  sequences  in  the  base 
were  selected  for  classification  studies.  The  parameters  for 
these  image  sequences  are  as  given  in  the  table  of  Figure  1.  This 
imagery  was  processed  as  follows: 

(1)  Windowing  to  eliminate  artificial  background  in  the 
original  data.  This  results  in  the  smaller  image  sizes  shown 
in  the  table. 

(2)  Optimum  global  binary  thresholding, basical ly  as  de¬ 
scribed  in  the  last  Final  Report.  This  scheme  assumes 
a  bimodal  Gaussian  distribution  for  the  image,  picking 
a  threshold  which  maximizes  the  between-cl ass  variance. 

The  technique  was  modified  to  disallow  thresholds  not 
contained  in  the  "saddle"  of  the  bimodal  distribution. 

The  threshold  ranges  obtained  by  this  automatic  routine 
are  given  in  the  table  of  Figure  1. 

(3)  Median  filtering  to  reduce  noise. 

Figure  2  contains  representative  pictures  of  the  eleven  groups 
after  preprocessing. 

The  images  thus  prepared  were  used  in  classification  studies 
utilizing  invariant  moment  (IM)  and  normalized  Fourier  descriptor 
(NFD)  features.  It  became  evident  that  the  three  groups  at  a 
distance  of  15  km.  (groups  3,  6,  11)  were,  for  practical  purposes, 
of  no  use 
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and  were  omitted  from  these  studies.  The  eight  remaining  groups 
were  represented  by  six  pitch/yaw  combinations  for  classifica¬ 
tion  purposes,  since  groups  one  and  four  differ  only  in  distance 
from  groups  two  and  five  respecti vely,  and  since  the  classifica¬ 
tion  is  invariant  to  size  change. 

The  groups  depicting  symmetric  views  (groups  4,  9,  10)  are, 
apparantlv,  slightly  nonsymmetri c .  It  may  be  that  a  bias  exists  in  the 
imaging  system  and  should  be  compensated  for  in  some  way.  Our 
experiments  did  not  explicitly  include  such  compensation. 

Synthesis  of  Imagery 

As  explained  above,  eleven  groups  of  twenty  images  were 
selected  from  the  NRL  data,  representing  all  the  distinct  pitch/ 
yaw/distance  combinations.  Neglecting  distance,  these  eleven 
groups  are  characterized  by  six  pitch/yaw  combinations.  For 
classification  purposes,  it  was  desired  to  obtain  an  "ideal" 
representation  for  each  of  the  pitch/yaw  groups.  Accordingly, 
as  an  approximation  to  these  images,  one  representative  image 
was  constructed  for  each  orientation  as  an  optical  (and  subjective) 
best  fit  to  the  data.  Pictures  of  the  synthesized  imagery  ap¬ 
pears  in  Figure  3. 

Unlike  the  real  data,  the  synthesized  Images  from  groups  4, 

9,  and  10  are  perfectly  symmetric,  since  they  are  intended  to 
depict  pictures  of  symmetric  objects.  The  features  from  these 
six  Images  were  used  to  characterize  the  six  pitch/yaw  combina¬ 
tions  for  classification  experiments. 
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Figure  3:  Symmetric  drone  images,  group  numbers  in  upper 
right  corner. 
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Feature  Extraction 


Feature  extraction  reduces  the  total  information  contained 
in  an  image  to  a  small  number  of  parameters.  These  parameters 
are  insensitive  to  change  in  size,  rotation  and  translation,  and 
are  easily  incorporated  into  a  classification  scheme.  Normalized 
Fourier  descriptors  (NFD's)  and  invariant  moments  (IM's)  were  the 
features  used  to  classify  the  NRL  data.  The  algorithms  are 
briefly  described  below. 


Normalized  Fourier  Descriptors 


The  Fourier  analysis,  and  the  associated  preprocessing,  are 
as  follows: 

(1)  Preprocessing  consists  of  threshold  and  median  filter. 

(2)  Following  the  median  filter,  a  chain  code  representation 
of  the  object  outline  is  constructed.  The  chain  code 
representation  constitutes  an  intricate  piece  of  soft¬ 
ware,  because  it  must  be  capable  of  handl i ng  a  number 

of  exceptions,  e.g.,  isolated  pixels,  object  extensions 
only  one  pixel  wide,  objects  with  embedded  "holes",  etc. 
The  ideal  would  be  the  usage  of  preprocessing  algorithms 
that  would  eliminate  such  anomalies.  However,  no  com¬ 
bination  of  preprocessing  steps  has  been  found  to 
entirely  eliminate  s  ~h  anomalies  under  the  variety  of 
of  combinations  of  r.oise  and  atmospheric  turbulence. 

An  alternative,  which  has  been  used  for  the  purposes  of 
studying  the  Fourier  descriptor  behavior,  has  been  to  use 
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the  Interactive  features  of  our  image  display  system  to 
"prune"  such  anomalies  by  manual  action.  Approximately 
seven  of  the  160  processed  NRL  Images  required  such 
acti on. 

(3)  Following  the  chain-coding,  the  chain-coded  elements 
are  interpolated  to  create  an  exact  power  of  two  number 
of  code  elements  for  Fourier  analysis.  The  FFT  is 
then  used  and  normalized  descriptors  are  calculated. 

(4)  Variation  of  the  number  of  retained  Fourier  coefficients 
from  N=16  to  N=512  was  carried  out.  From  these  varia¬ 
tions  it  was  established  that  a  small  number  of  co¬ 
efficients,  e.g.  16_<  N<64,  are  sufficient  to  produce 

a  good  representation  of  the  object  outline. 

The  NFD's  normally  constitute  an  excellent  characterization 
of  the  object.  However,  in  the  presence  of  noise,  there  is  a 
possibility  of  an  irremediable  error  in  the  normalization  calcu¬ 
lation'  Noise  in  the  original  image  causes  perturbations  in 
the  corresponding  Fourier  coefficients.  In  particular,  when 
the  coefficients  are  of  small  magnitude,  the  noise  can  cause  a 
sign  change  in  the  real  part  of  those  coefficients.  Since  the 
normalization  criterion  is  to  maximize: 

?  Ret F(k) > I F(k)  I 
k 

normalization  errors  can  result.  For  group  9,  this  feature  error 
was  significant  and  caused  significant  classification  error.  The 
literature  contains  no  mention  of  this  problem. 
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Invariant  Moments 


The  analysis  by  invariant  moments,  and  the  associated  pre¬ 
processing  are  as  follows: 

(1)  Preprocessing  consists  of  threshold  and  median  filter 
as  in  the  NFD  calculations. 

(2)  Following  the  median  filter,  the  edges  of  the  object 
are  determined  as  in  the  bounda ry- traci ng  algorithm 
associated  with  the  NFD's  (see  Final  Report,  1982). 

(3)  The  seven  IM's  of  the  object  and  the  object  boundaries 
are  computed,  giving  fourteen  numbers  invariant  under 
size  change,  rotation  and  translation.  The  inclusion 
of  the  seven  moments  of  the  edges  is  to  insure  good 
description  of  objects  with  prominent  edges,  but  which 
occupy  little  area. 

Classification  Libraries 

A  classification  library  is  a  set  of  features  computed  from 
objects  known  a-priori.  Given  a  piece  of  data,  the  solution  of 
the  classification  problem  is  to  minimize  a  criterion  function 
with  respect  to  the  elements  of  the  library.  The  nearest- 
neighbor  criterion  gives  a  simple  and  useful  example  of  such  a 
functi on. 

The  objects  in  the  NRL  data  are  the  noisy  two  dimensional 
images  of  the  drone  at  various  orientations.  Since  no  perfect  images 
were  available  to  construct  a  classification  library,  two 
alternatives  were  pursued.  The  first  was  to  use  the  features 


of  the  six  synthesized  drone  images  as  a  library.  These  features 
correspond  to  the  noiseless  imaging  case  and  are  referred  to  as 
the  ideal  library. 

The  second  alternative  was  to  construct  an  "average-data" 
library  which  likewise  contained  six  elements  but  now  each 
element  is  the  group  average  of  the  features.  It  was  reasoned 
that  such  a  library  should  be  more  accurate  in  the  presence  of 
noise  and  would  partially  compensate  for  any  bias  in  the  imaging 
system.  Both  libraries  were  used  in  the  classification  experi¬ 
ments  . 

C 1  ass i fi ca tion 


Classification  based  on  a  nearest  neighbor  criterion  is  the 
simplest  and  was  considered  first.  Let  N  be  the  dimension  of  the 
features.  For  IM's,  N*14,  while  for  NFD's,  N=512.  Denote  the 
feature  vector  of  a  piece  of  data  by 


D  •  <Vi-i 


and  the  feature  of  the  Mth  class  of  a  library  by 


The  classification  of  the  data  is  accomplished  by  determining 
that  value  of  j  for  which: 


|CJ-D[  Ip^.!  jCK-Dj  |p  for  all  k.  The  data  is  then 


"closest  to"  and  is  considered  to  be  a  member  of  class  j 

N 


1  P 


is  the  £p  norm  given  by 


11*11 


"real,  &  • 


This  norm  with  p  =  1 ,2  was  used  to  classify  the  NRL  data. 
Classification  by  NFD's  is  described  first. 
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Nearest  Neighbor  Classification  by  NFD's 

Nearest  Neighbor  Classification  by  NFD's  was  carried  out  with  16,  32  and 
64  retained  coefficients.  The  results  are  presented  in  the 
"confusion  matrices"  of  Figures  4-9.  The  matrices  demonstrate 
that  the  classification  is  very  insensitive  to  change  in  the 
number  of  retained  coefficients  for  16<N<64.  In  fact,  there  is 
almost  no  advantage  in  choosing  more  than  16  coefficients.  It  is 
therefore  clear  that  N  = 1 6  is  the  best  choice  for  economy  of 
representati on . 

2 

For  the  NRL  data,  the  use  of  the  l  norm  resulted  in  fewer 
errors  as  can  be  seen  from  the  following  table  derived  from 
Figures  4-9: 


Norm 

Library 

%Error 

.e1 

averaged  data 

32 

i  deal 

29 

31  average  error 

CM 

averaged  data 

31 

CM 

i  deal 

24 

28  average  error 

Table  1:  Average  Error  for  NFD  Classification. 
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Figure  6:  Nearest  Neighbor  NFD  Classification  with  averaged- 
data  library  and  64  retained  coefficients. 

Upper  =  O  norm,  lower  *  l 2  norm.  Heavy  outline 
represents  correct  decision. 


Figure  7:  Nearest  Neighbor  NFD  Classification  using  ideal 

library  and  16  retained  coefficients.  Upper  *  .C1 
norm,  lower  *  £2  norm.  Heavy  outline  represents 
correct  decision. 
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Figure  8:  Nearest  Neighbor  NFD  Classification  with  ideal 

library  and  32  retained  coefficients.  Upper  =  i 
norm,  lower  =  £2  norm.  Heavy  outline  represents 
correct  decision. 


In  compiling  this  table,  group  9  was  omitted  to  insure  that  norm¬ 
alization  errors  did  not  affect  the  choice  of  norm.  Despite 
the  30%  or  so  error  listed  in  Table  1,  the  NFD  Classification 
is  an  excellent  one.  Examination  of  the  confusion  matrices  of 
Figures  4-9  and  examination  of  the  original  images  reveals  that, 
while  the  classified  group  may  not  always  be  the  correct  one,  the 
selected  group  always  bears  a  strong  resemblance  to  the  correct  one. 
That  is,  perturbations  in  the  NFD  coefficients  correspond  in  a  natural 
way  to  perturbations  in  shape. 

It  is  interesting  to  note  that  use  of  the  ideal  library  re¬ 
sulted  in  a  lower  percentage  error  than  when  the  averaged-data 
library  was  used.  Since  a  noise-free  library  may  be  the  only  one 
available  in  practice,  this  fact  makes  the  NFD  Classification 
attracti ve . 

In  summary,  the  confusion  matrices  of  Figures  4-9  indicate 
the  use  of  16  retained  coefficients  with  the  norm  while 
Table  1  supports  the  use  of  an  ideal  noise-free  library.  The  Nearest 
Neighbor  NFO  Classification  gives  excellent  performance. 

Nearest  Neighbor  Classification  by  IM's 

The  results  of  Nearest  Neighbor  Classification  by  IM  are 
presented  in  the  confusion  matrices  of  Figures  10  -  11.  For  the 
NRL  data,  the  use  of  the  norm  resulted  in  slightly  fewer  errors,  as 
examination  of  the  following  table  reveals: 
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Norm 

Li  brary 

%  Error 

averaged  data 

41 

i  deal 

54 

48  average  error 

CM 

averaged  data 

53 

*2 

1  deal 

59 

56  average 

Table  2: 

Average  Error  for  Nearest  Neighbor 

IM  Classification 

However,  use  of  the  norm  improves  the  situation  only 
slightly  and  the  performance  of  the  classification  is  disap¬ 
pointing  overall.  The  average  error  is  about  twice  as  much  as 
in  the  NFD  classification  with  16  retained  coefficients.  It  is 
possible  that  the  low  performance  of  an  INI  classifier  ^s  due  to 
large  correlations  between  the  moments,  in  which  case  the  near¬ 
est-neighbor  criterion  is  suboptimum.  This  question  is  taken 
up  later. 

Comparison  of  Nearest  Neighbor  Classification  Schemes 

The  performance  of  the  classification  by  NFD's  was  consider¬ 
ably  higher  than  that  by  IM's,  the  NFD  average  classification 
error  being  about  1/2  that  of  the  IM.  Furthermore,  the  NFD 
classification  was  less  affected  by  imaging  system  bias  as  can 
be  seen  from  Tables  1  and  2.  Table  1  indicates  the  NFD  error 
is  lower  with  the  ideal  (unbiased)  library  whereas  Table  2 
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Figure  10:  Nearest  Neighbor  IM  Classification  with  average 
data  library.  Upper  =  l 1  norm,  lower  *  £2  norm 
Heavy  outline  represents  correct  decision. 


Figure  1 1  : 


Nearest  Neighbor  ^  Classification  with  ideal 
library.  Upper  =  V  norm,  lower  =  l 2  norm. 
Heavy  outline  represents  correct  decision. 


indicates  IM  error  is  lower  with  the  a veraged-data  (biased)  li¬ 
brary.  These  facts  make  the  NFD's  the  feature  of  choice  over 
IM 1 s  for  NN  classification. 


Maximum  Likelihood  Classification  by  IM's 

In  an  attempt  to  improve  the  performance  of  the  IM  clas¬ 
sification,  a  criterion  was  chosen  to  reduce  the  effect  of 
moment  correlations.  Specifically,  the  approach  may  be  stated 
as  follows.  The  features  of  a  piece  of  data  are  known  to  con¬ 
tain  one  of  M  library  features  C1  plus  additive  noise.  The 
outcomes  are  assumed  equally  likely.  That  is,  the  ith 
hypothesis  is: 

H.:  d  =  C1  +  n  where  d.C^.n  are 
N  dimensional  vectors,  the  Ci  are  equally  likely  and  n,  the  noise, 
is  zero  mean,  and  Gaussian  with  covariance  matrix  Rn.  Given  a 
piece  of  data  D,  the  optimal  classification  qrouD  (P  minimizes 

(D  -  Ci)TRn_1(D  -  C1) 

over  all  elements  of  the  classification  library  [1],  If  the 
inverse  does  not  exist,  a  pseudoinverse  must  be  used.  Numeri¬ 
cally,  a  pseudoinverse  makes  sense  even  if  the  matrix  is  only 
i 1 1 -condi ti oned. 

The  global  covariance  matrix  of  the  IM  data  features  was 
computed  and  is  listed  in  Figure  12  (the  matrix  is  symmetric). 

This  matrix  was  computed  as  follows: 

(1)  The  group  mean  was  removed  from  each  of  the  eight 
data  groups  to  yield  a  zero  mean  vector  sequence 


(2) 


.  k, 
{n  } 


160 


k=l  (160  is  the  number  of  pieces  of  data). 

The  ijth  element  R.  .  of  the  covariance  matrix  was  calcu 

*  J 

lated  using  the  formula: 


„k  „k 
n.  nj 


It  can  be  seen  by  inspection  that  Rn  is  nearly  singular,  hence 
ill-conditioned.  Accordingly,  the  Moore-Penrose  pseudoinverse 
was  used  for  the  required  inversion.  The  pseudoinvers  is  listed 
in  Figure  13  and  the  maximum  likelihood  classification  based  on 
it  is  summarized  in  Figure  14. 

The  performance  is  only  trivially  better  than  the  Nearest  Neigh¬ 
bor  IM  Classification.  Table  3  presents  the  average  error. 

Li  brary _ %  Error 

average  36% 

ideal  57% 

46%  average  error 


Table  3:  Average  Error  of  Maximum  Likelihood  IM  Classification. 

This  46%  is  to  be  compared  with  the  48%  from  Table  2,  and  is 
still  far  above  the  28%  listed  in  Table  1.  Thus  no  IM  classifi¬ 
cation  is  nearly  as  good  as  the  Nearest  Neighbor  NFD  Classification 
with  16  retained  coefficients. 

ine  probable  reason  for  the  poor  results  of  Figure  14  and 
Table  3  is  that  the  noise  is  not  stationary.  Specifically,  the 
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features  of  some  groups  exhibited  much  more  variation  than  others. 
The  covariance  matrix  is  a  global  measure,  however,  and  cannot 
account  for  this  nons tati onari ty . 

Conclusions 

The  preprocessing,  feature  extraction  and  classification 
algorithms  have  been  applied  to  classify  eight  groups  of  twenty 
images  taken  from  the  NRL  data.  The  features  used  in  the  classi¬ 
fication  experiments  were  the  normalized  Fourier  descriptors 
(NFD's)  and  the  invariant  moments  (IM's).  Use  was  made  of  two 
classification  libraries:  the  ideal  library  comprised  of  features 
derived  from  synthesized  imagery  and  the  averaged-data  library 
consisting  of  the  group-averaged  features.  The  Nearest  Neighbor-NFD 
classification  with  16  retained  coefficients  was  the  best  overall 
performer  in  terms  of  accuracy  and  economy  of  representation. 

The  best  IM  classification  was  only  1/2  as  accurate  as  the  six¬ 
teen  coefficient  NFD,  in  an  average  sense. 

NFD's  are  normally  an  excellent  characterization  of  an  object. 
However,  they  are  subject  to  two  problems:  noise  in  the  original 
image  can  cause  normalization  errors,  and  the  algorithm  requires 
a  sophisticated  boundary- traci ng  subalgorithm.  Both  of  these 
topics  merit  more  research.  These  facts  notwithstanding,  the 
NFD'S  are  well-suited  to  classification  of  aircraft  in  turbulent 
atmosphere . 
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Appendix  -  Software  Listing 

(1)  Contr:  Boundary-tracing  routine. 

(2)  Interp:  Interpolates  a  sequence  to  a  power  of  two. 

(3)  FD:  Normalizes  the  Four  descriptors. 

(4)  MOMNT:  Calculates  the  seven  invariant  moments  of  an  image. 


IF  •'  LIME  *  GE  ,  MM  LIME  >  GO  TO 

rt  p  2  0  2  J  :  \  ••  m ;  a  3  ,A,  ;.-i  p 

COMT [HUE 
run  rn  iai 


E  THE  y 0 R K F I L  E 

C  A  L  L  POST*'  W  F % 

KDUHT  KOUMT  —  1 

_t%  THE  LIME  ADJACIEHCY  SCRATCHFIi 


OF  EM  ( UM  JT~  i  i  ?  NAME58 '  LAG  ’  ?  TYPE  - ' ; 
rj  n  too  is  -  1  ..  q  jj  j\i  T 

VJP I TE  ( 1 1 '  K  >  <  I A  v  J )  ?  J  - 1 7  4 ) 


CCMTIMUE 

STERMIMS  WHICH  LINE  SE3MEMTS  ARE  APJACIEMT*  STGRIr 
ESPECTIUE  KTH  LIME  ABJACIEHCIES  IM  THE  KTH  RECORD 

TJ  0  4  0  0  M  -  2  ?  K  0  U  M  T 
■!  ..1 

-  >  ■  W  W  ■  •/  .U  j  t 

T •:%  n  !  *'  I  '»  :r  A 

•«.  .  I  ««•  >  to/  .*  W 

COr!  T I  ;•!'.!£ 

T  3  T  ©  T  ::=  M  ft  X  0  (  1  f  <  M  —  1  2  )  > 

DO  -101  K  *  1ST  ART  ?  M-l 

I F ■ ABS <IL< K > - IL ( N >  > . GT . 1 . OR . IL  < K  >  *  EQ . I L < H ) >  30 
IF  < 1X2 < K > . LT , ( IXi ( N)-l > , OR  *  1X1 < K ) . GT , < 1X2 <  H >  +  i: 
READ  ( 1 1 '  K  >  <  I A  <  J )  y  J=1  ?  <1 ) 

MUM  -  MUM  +  1 

.00  ■I 03  J  =  1?4 

IF  <  ZA<  J  >  » ME  *  0  >  GO  TO  -103 

I.AO J  <  HUM  >  ~  K 
30  TO  />,0'5 
CONTINUE 
C OMTIMUE 

UR  I TE < 1 1 ' K >  < I A { J  >  ?  j- 1 ?  4) 

WF:  I  TE  <  1 1 '  H )  <  I  AD  J  ( J ) .  J-- 1?4) 

COMTIMUE 
0 OMTIMUE 

TART  THE  CONTOUR  TRACING  ALGORITHM 


-JMM  -s  1 


■ETERM I  ME  WHICH  SEGMENT*  ARE  AD  AC  TENT  TO  THF  FIRST  D 


iz’JrHzr!  l  s  ;  iRL  ADLALILNi 


READ  (  i  1 ' 1  >  \  I A  ( J  >  ?  J  - 1  -  -1  > 


ETERMIHE  THE  DIRECTION  OF  THE  FIRST  LIMF  SEGMFHT 


A  T 


Tff  KJHT 


•=  p«j 

M 

+  i 

I  HUE 

HE  E 

NT 

IRE  CONTOUR 

HAS  BEEN 

TRACED 

A  »  c 

( A 

( 1 } -A ( MUM 1 ) > 

A  03 

( B 

i  >  -  » .*  NUM1  >  > 

LOOK 

sr 

OR  SOMETHING 

ELSE  TO 

DO 

T  #  j  £ 

.  1 

.  AND . 3T , LE . 1 

♦  ;'ir 

IT:  _  MI  IW 

,  GT ,  ■! )  GO 

{ 1 1 / 

II 

N> (IA< J>  t J=1 

,  4  ) 

I  IN' 

TH 

OLD  FLAG 

rn  4  9P 


c  •:  1 1  m  > 


rrt 


i.: 

[j 

r 

0 

B 

H 


LOOK  FOR  THE  CLOSEST  NEW  LINE  SEGMENT  WITH  RESPECTIVE  DIRECTION 

DO  603  K  =  It  A 
TINT  «  I A  ( K ) 

IF  \  .TINT  ,  EO  .  0 )  GO  TO  603 
IF ( C ( I INT  >  ♦  NE . 0  >  GO  TO  603 
ID15;:ABS  < I  NT  ( A  ( M  U  M 1 ) )  - 1 X 1  <  1 1 H  T )  > 

ID2=ABS ( INT (A < HUM! ) ) -1X2 { I INT ) ) 

IF •  I.D1  .OS.  ID.  AND .  ID2 . GE .  ID )  GO  TO  603 
ILL  --  IL<  IIMT> 

I IN  «  I INT 

IF < ID l . GT , ID2  >  GO  TO  605 
IXX1  -  1X1 < I. INT ) 

I XX  2  1X2  <  1 1  NT ) 

r  n  ~  i  ri  i 

INC  a  1 
GO  TO  603 
CONTINUE 
I XXI  a  1X2(1 INT) 

IXX2  =  IXi < I. I  NT) 

ID  ■•=  ID2 
INC  a  -1 
CONTINUE 

IF  NO  NEW  SEGMENTS  TO  BE  FOUND ?  BRANCH  TO  TRACES ACK  SECTION 

IF  ( ID  ♦  EO  >  1000)  GO  TO  7J?9 
NUML  —  NUN!  -  1 


r  r  ( ;‘.53  i;.jt  {  A  .  ;:•*?  .'ML  >  gr  .  £  >  nr.  rn  ;3i 

-  TihiiMarr'  r  '-I*81  "-liY 


A 


rn 


GO  TO  SOI 


r  p :  r  ;>•;  \  .  )  0  q  t  q  4  7  » 

00  TO  6 SO 
;7~  CONTINUE 

•3<rj.  CONTI  HUE 

I?  NO  NEW  3 (EON ENTS  -  00  r0  EN.D  SECTION 
r  p  j.jsry  .  ,M£  4  1  )  (30  TO  9 1 1 

<\{  IJ  l‘j  .«  j\l  1 1  j  ^ 

TRACE  BACK'  SECTION 

.31  CONTINUE 

•  HUM  )  a  A  <  NUNL  > 

3- HUM)  =*  3  ( NUML  > 

NUM1  -  HUM 
NUN  =  HUM  r  1 
NUML  »  NUML  -  1 
10  =  1000 

SEE  IE  THERE'S  NON  AN  A 3 J AC I ENT  SEGMENT  TO  THE  U 

DO  £90  K  »  1 ? KQUNT 

I F  <C( K ) . HE . 0 )  GO  TO  690 

T.  D  .1  *  ABS  <  I  NT  <  A  <  HUM  1 ) )  -  IX 1  <  K  )  ) 

132  ■=  ABS  <  I  NT  ( A  ( HUM!  >  )  -1X2  ( K  >  ) 

I 03  «  ABS ( I NT (8 < NUM 1)>- I L ( K > > 

IF  (  <  I D 1 .  OT ,  1 ,  AMD ,  ID2  ,  GT ,  1 }  ♦  OR  ♦  I  Do .  GT ,  1  >  GG 
IF  '  ID! , GE , ID . AND . ID2 , GE * ID >  GO  TO  690 

DETERMINE  CLOSEST  NEW  ADJACIENT  SEGMENT  WITH  D.l 

.TIN  «  K 
ILL  a  ILCK) 

IF<ID1.GT*ID2>  GO  TO  691 
IXX1  a  IX 1 <K> 

IXX2  =  1X2 (K) 

INC  L 

GO  rr'  -ac. 

;  a  •:  rj  n  m  t  t  m  j  j  jr 

I XX I  -  1X2 (K) 

IXX2  --  1X1(10 
INC  -1 
ID  a  ID2 
:  9 0  CONTINUE 

NEED  TO  TEACE3ACK  MORE  ? 

IF  <  ID  »EQ .  1 0 0 0  >  0 0  T 0  6 3 1 

nr.  rr;  4 r.  ; 


SECTION 

rs.viT'ra  rnr 


i  r  <  I'lJt’c  (  A  <  NUsi'l  >’  -Vi  \  )  )  lijt  ♦'!  *  UK  ^  i-i.b;s 

M  =  K 

‘30  TO  6? s 

CONTI HUE 

TP AC SB A CK  TO  BEGINNING 
CONTINUE 

2r(ii*EQ.  I  >  00  TO  a?? 

DO  697  K  =  M?2?-l 
A ( MUM >  =  A ( K > 

3 ( NUN >  -  9(K) 

MUM1  -  MUM 
MUM  •-  MUM  +  1 
CONTINUE 
CONTINUE 

DELETE  SCEATCHFILE 
CLOSE < UNIT-11 ) 

OPEN  < UNIT  -11 9  NAME  - 'CONT , CUT ' > 
WRITE! 117511)  NUM1 
DO  307  K  ~  IfNUMl 
WR I TE  ( 1 1 7 51 2)  A ( K > ? B ( K ) 

CONTINUE 
FORMAT (2X7 13) 

:  FORMAT (2Xr2F6.0> 

CLOSE (UNIT-11 > 


u'J  IU  ar'J 


SET  UR  THINGS  FOR  THE  IIS 

ACCEPT  977 » ANS 

IF ( ANS . NE  *  '  Y '  >  GO  TO  999 

TYPE  500 

ACCEPT  977 ? ANS 

FORMAT (A2> 

IF (ANS, EG, 'R ' )  GO  TO  91 3 
CALL  GRINT(IXl) 

»•  j  g  |  jv--  ^ 

I XI (K)  =  -1 
CONTINUE 
NCHN  »  -32763 
TYRE  501 
ACCEPT  *7  NMA3K 
GO  TO  914 
l  NMASK  =  255 

TYPE  501 
ACCEPT  7  7 NCHN 
\  CONTINUE 

SCALE  THE  COORDINATE  FOR  THE  IIS 

DO  933  K  -  1 7 MU Ml 
A ( K )  -  A ( K )  -  M MS AMP/ 2  v  255 
B(K)  -  B ( K )  -  MMS AMP/2  -f-  255 
:  CONTINUE 

WRITE  THE  CONTOUR  TO  IIS 


X*f  (J  (2 


I  ,  ;  M;  HA  1  _  •! 


t  r.j  r  {  u*;  \  \ 


£  0  c 

FOR  Mi1)  T  <’  A2  ,* 

93? 

CONTINUE 

9  4 1 

CONTINUE 

C  r 

L uSE  f M E  0 0 N 7 0 U R 

IXXi  ~  INT<A<HL<M1 )  > 

I VI  -  [HTiSiNUMl ) > 

IXX2  -  I NT  ( A  •:  i  >  > 

IV. 3  ~  INTO <  1  >  ) 

CALL  DVECT  <  FC?  ?  IXXI ? IYi  » 

IXX2  ? I Y2 ?  NCHN ?  NMASr 

occ 

CONTINUE 

C  i-f 

LL  DONE 

- 

STOP 

•*%  i\  (\ 

FORMAT <  1 1 X  ?  ' DO  YOU  WANT 

GRAPHICS  OR  REPRESS- 

FORMAT < 1 IX?' WHAT  COLOR  ? 

( 1--REB?  2-GREEN?  > 

1  •/*  «“ 

FORMA  Til  IX  r  '<30  ?  '  ?  9  ) 

5 1 0 

FORMA T  f 1 1 X  ? ' DO  YOU  WANT 

DISPLAY  ?  '?$) 

36 


I 


I 


SUBROUTINE  IHTERP • CHNR * CHHC > 
DIMENSION  CHNR  <  1024  > > CHNC  <  1021 
DIMENSION  NAME ( 15) 

DATA  NAME/ ' CO ' ? ' NT  ' r  ' UR ' ? ' . 0  '  r 
SQ2  »  SORT (2.) 


)  ? RTEMP <  102  1 )  7 CTEMP <  102-1) 
' UT 9*0/ 


***********************************************  ******************* 


*****  OPEN  OLD  FILE  AND  READ  THE  CHAIN  CODE  VALUES . 
*****  NTOT-  TOTAL  I  OF  POINTS 

*****  CHNR CI)~  REAL  PORTION  OF  ITH  CHAIN  CODE 

*****  CHNC ( I )  -  IMAGINARY  PORTION  OF  ITH  CHAIN  CODE 


TYPE  11 


.  11  FORMAT ( 1.1  X»  'ENTERING  INTERP  '  ) 

I  TYPE  12 

1  12  FORMAT  ( .1 1 X  ?  '  ENTER  THE  INPUT  VERSION  NUMBER  ..  '?$> 

ACCEPT  13 » NAME (7) 
l  12  FORMAT (A2> 

i  OPEN ( UNIT=1 1 , NAME-NAME 7  TYPE- ' OLD ' > 

READ  < 1 1 r 1 ) * NTOT 
DO  10  I”1 » NTOT 

!  READ ( 1 1 7 2 ) RTEMP ( I ) » CTEMP < I ) 

.10  CONTINUE 

1  FORMAT <2X» 13) 

I  2  FORMAT < 2.X, 2F6.0) 

•  CLOSE  <  UN I T~ 1 1 ) 


■  C 

,  C  ***** 
C 


t 


1 


* 

1 

1 


F 

4,c  ***** 


c 


K -1 

DETERMINE  CONTOUR  LENGTH  AND  DISTANCE  BETWEEN  POINTS 


CONLEN-O 
DO  25  1=2  r NTOT 
RING  =  1 

Rl-RTEMP ( I ) -RTEMP ( 1-1 ) 

Cl  =CTEMP  ( I )  -CTEMP  ( I-.1 ) 

IF ( R1 . NE . 0 . AND . Cl . NE . 0 )  RINC  ~  S02 
C  0  N  L  E  N  =  0  0  N  L  E  N + R I N  C 
CONTINUE 
RINC  =  1 

Rl-RTEMP ( 1 ) -RTEMP ( NTOT ) 

Cl =CTEMP ( 1 ) -CTEMP ( NTOT ) 

IFCR.l  . HE, 0. AND, Cl  ,NE.O)  RINC  =  S02 
CONLEN=CONLEN+RINC 
DIV-CONLEN/512 ♦ 

STARTING  PROCEDURE 

PLEN  »  1 

RL-RTEMP ( 1 ) -RTEMP ( NTOT ) 

CL-CTEMP ( 1 ) -CTEMP ( NTOT ) 

I F  <  RL ♦ NE  *  0  *  AND ♦  CL » ME . 0  >  PLEN  =  S02 
RMULT =K*DIV 

IF (R'MULT .  GT .PLEN) GO  TO  <10 

FACT--RMULT/PLEN 

CHNR ( K ) = RTEMP ( NTOT ) +FACT*RL 

CHNC  ( K )  ~CTEMF'  ( NTOT )  +FACT*CL 

K-K+l 

C-0  TO  30 

CONTINUE 


-J. 


I SAVE* 1 

R.l  ~RTEMP  ( 1 )  -RTEHP  < NTOT ) 

Cl -CTEMP ( 1 > -CTEMP <  NTOT ) 

PPLEN  »  1 

I F  ( R1 * NE . 0 . AND .Cl. NE . 0 )  PPLEN  =  SQ 
CONTINUE 

RL  =  RTEMP ( ISAVE+1 >-RTEMP ( I SAVE) 
CL=CTEMP ( ISAVE+1 > -CTEMP ( ISAVE > 

PLEN  =  1 

I F  <  RL ♦ NE . 0  * AND ♦  CL ♦ NE . 0  >  PLEN  *  SQ2 

CONTINUE 

RMULT  -  K'fcDIV 

IF { K ,  GT , 512 )  GO  TO  90 

I F • RMULT . GT . < PPLEN+PLEN > )  GO  TO  S5 

FACT  *  (RMULT  -  PPLEN) /PLEN 

CHNR ( K )  *  RTEMP< ISAVE)  +  FACTKRL 

CHNC(K)  *  CTEMP (ISAVE)  +  FACTKCL 

K=K+1 

GO  TO  50 

ISAVE -ISAVE+1 

PPLEN  »  PPLEN  +  PLEN 

GO  TO  60 

CONTINUE 

TYPE  14 

FORMAT <  1 IX  r 'LEAVING  INTERP ' ) 

RETURN 

END 


s UBRD U T I H E  FD ( CHNR  ? CH M C?HS> 


r> 


CALCULATES  THE  FOURIER  DESCRIPTOR  OF  A  CONTOUR 

REAL  CHMR  <  5 1 2  >  ?  CHK'C  C  51 2  > 

REAL  MXR ? M I R ? MXC ? M I C 
REAL  MR t MC 

REAL  CHMTR ( 512 )  ? CHNTC < 512 >  mMAGI  r  MAC-2. -MAC- 
DOUBLE  PRECISION  ARC  ?  ARCT  ?  SR  ?  P I  ?  P 1 2  ?  !J  ? V 
INTEGER  FCB ( 6 >  ? BUFF ( 236 ) 

DATA  PI/3 . 1 T 1 59265338979323S/ 

p  T  ">  —  2>kP  I 

D  A  T  A  9  U  F  F  / 2 5 6 * - 1  / 

DATA  N/512/ 

rypc  ;v .  '  ENTERING  FD  ' 

r . •. r  FFT 

CALL  FFT2  <  CHMR  ?  CHNC  r  N - 1) 

SEE  IF  THE  CONTOUR  IS  TRACED  IN  THE  CLOCKWISE  DIRECTION 


MAG2 


SORT  <  CHNR < 2 ) ** 2+ CHNC < 2  > > 
one t / CHNR <  N >  $$2+ CHNC < N ?  2 > 


IF  NOT?  REVERSE  THE  SEQUENCE 

IF<MAG1.GE.MAG2>  GO  TO  321 
DO  311  K  ~  2?N 
CHNTR(K)  -  CHNR(K) 

CHMTC < K )  -  CHNC(K) 

3 11  CONTINUE 

DO  312  K  ~  2  r  N 
ruwDfk;)  -  CHNTR  (N+2-K  > 

CHNCCK)  -  CHNTC ( N+2-K ) 

312  CONTINUE 
MAGI  -  MAG2 

321  CONTINUE 

MAGI  IS  THE  MAGNITUDE  OF  F  C 1  > .-  U  IS  THE  PHASE  OF  F(1 

U  "  ATAN( CHNC ( 2 ) /CHNR ( 2 } ) 

IF ( CHNR  (2 )  *LT . 0 )  U  =  U  —  F'l 


NORMALIZATION 
NORMALIZE  POSITION 


CHNR  (  1  \ 
CHNC  •’  1 > 


A 

o 


MAG2  ■“  0 
DO  200  K  ■»  2  ?  N 

DIVIDE  BY  THE  MAGNITUDE  OF  THE  LARGEST  COEFFICIENT 


CHNR  < K /  -  CHNR < K ? /MAG 1 

CHNC ( K )  -  CHNC \K> /MAGI 
MAG  '  SORT { CHNR  <K>  :}::J:2  f  CHNC  ( K )  :{c:S2  ■ 


41 


IF ( K . EG . 2 . OR . MAG .  LE . MAG2 >  GO  TO  200 

HA02  —  M G 
-  |; 

CONTINUE 


C  F < K2 >  IS  THE  COEFFICIENT  OF  SECOND  LARGEST  MAGNITUDE 
C 

C  V  IS  THE  PHASE  0F'F(K2) 

C 

V  ~  ATAN ( CHNC C  K2  > /CHNR <  K2 ) > 

IF- CHNR(K2> .LT.O)  V  ~  V  -  PI 

C 

C  SHIFT  l\2  TO  BE  IN  (-N/2+l»N/2> 

C 

K2  K2  -  1 

IF ( K2 . GT . N/2 >  K2  -  K2  -  N 
C 

C  CALCULATE  THE  NORMALIZATION  MULTIPLICITY  OF  F<K2> 

C 

MK2  “  ABS i K2-1 ) 

C 

C  PERFORM  THE  FIRST  OF  THE  MK2  NORMALIZATIONS 

C 

DO  370  K  ~  2 ?N 
K1  -  K  -  1 

IF  ( K1 . GT . N/2 )  K1  •=  K1  -  N 

SR  ^  <<K1-K2>*U  +  < 1-K1 )  *V  )  /FLOAT ( K2-1  > 

CHNTR(l)  *  CHNR(K)*DCOS(SR)  -  CHNC\K>*DSIN<SR> 

CHNC  <  K  >  ■=  CHNR  <  K )  KDS I N  <  SR )  +  CHNC  (K  >  *DCOS  < SR  > 

CHNR ( K >  "  CHNTR C 1 ) 

CHNTR(K)  =  CHNR ( K  > 

CHNTC (K>  ~  CHNC(K) 

375  CONTINUE 

C 

C  CALULATE  THE  CORRESPONDING  ABIGUITY  RESOLVING  CRITERION 

C 

ARC  ~  0 

DO  582  K  =  2  ? N 

ARC  =  ARC  +  CHNR  <  K ) YABS  <  CHNR  <  K )  > 

5S3  CONTINUE 

C 

C  PERFORM  THE  REST  OF  THE  NORMALIZATIONS  WITH  THEIR  RESPECTIVE  ARCS 
C 

IF <MK2 .EG . 1 )  GO  TO  901 
DO  396  NUM  -  2  7 MK2 
ARCT  =  0 
DC  33-1  K  -  2  7  N 
K1  -  K  -  1 

IF ( K1 . GT . N/2 )  K1  -  K1  -  N 
SR  »  <  K 1  - 1 )  :KP  1 2/FLOAT  <  MK2  > 

SR-  -  SR 

CHNTT  -  CHNR<K>*DCOS(SR>  -  CHNC(K)*DSIN(SR> 

CHMC(K)  -  CHNR  < K ) YDS I N < SR )  +  CHNC(K)*DCOS(SR> 

CHNR  <  K )  -  CHNTT 

ARCT  -  ARCT  +  CHNR  < K ) #ABS ( CHNR { K ) > 

C 

C  PICK  THE  NORMALIZATION  WITH  THE  LARGEST  ARC 
C 

OS'!  CONTINUE 

IF ( ARCT. LE. ARC)  GO  TO  386 
DO  333  K  ^  2  r  N 
CHNTR ( K  )  ■=  CHNR(K) 

CHNTC (K)  -  CHNC ( K > 

rnMT  r  mi  tc 


i 


<  II  k  W  I 


'  *'  *  •  »l 

586  CONTINUE 

DO  902  K  -  2yN 
CHNR  <  K )  *  CHNTR ( K > 

CHNC(K)  »  CHMTC(K) 

902  CONTINUE 

C  FILTERING  -  IDEAL  LOPASS 

c 

"*901  IF (  (NS+1 )  . GT.N)  GO  TO  555 

NC  -  NS/2  +  1 
ND  =  N  -  NS/2 
DO  1  K  ~  NC  t ND 
CHNR(K)  »  0 
.CHMC(K)  ®  0 

1  CONTINUE 

-*k 

C  TAKE  INVERSE  FFT  OF  THE  FOURIER  DESCRIPTOR 

C 

555  CONTINUE 

OPEN  ( UNIT-- 11? NAME  =  '  FD .  OUT  '  t FORM® '  UNFQRMATTf 

DO  391  L  *  1 ?  312 

WRITE (11)  CHNR ( L ) ? CHNC ( L > 

G91  CONTINUE 

CLOSE (UNIT-11) 

CALL  FFT2 ( CHNR t CHNC rNf+l) 

DO  2  K  ”  lrN 

CHNR < K >  ~  CHNR <K> /FLOAT <N) 

CHNC<K)  *  CHNC (K) /FLOAT (N) 

2  CONTINUE 

TYPE  502 
ACCEPT  503  y ANS 
IF ( ANS .NE .  '  Y  '  >  GO  TO  333 
C 

C  SET  UP  THINGS  FOR  THE  IIS 
N  =  512 


TYPE  500 
ACCEPT  IOOjANS 
100  FORMAT (A2> 

IF(ANS.EQ. 'R'  >  GO  TO  201 
CALL  GRINT (BUFF ) 

DO  121  K  =  ly256 
BUFF(K)  ■-  -1 
121  CONTINUE 

NCHN  ■  -32763 
TYPE  501 
ACCEPT  M,  NMASK 


GO  TO  202 
201  NMASK  *  255 

TYPE  501 
ACCEPT  M ? NCHN 
20  2  CONTINUE 

SCALE  COORDINATES  FOR  IIS  DISPLAY 


}  -  -1000 
l  =  1000 
:  =  -loco 
:  ■  iooo 

601  K  ■»  1 1 M 
5  =  MAX ( MXR 7 CHNR  < 
i  =  MIN (MIR  7 CHNR < 
:  ■=  MAX  (MXCy  CHNC  < 
:  -  MIN (MIC ? CHNC ‘ 


MC  ~  MXC  -  MIC 
MIR  ■=  MIN  (MIR?  MIC) 

SC  =  4  12/MAX  <  MR  ?  MC  > 

DO  602  K  -  I ?  N 

CHMR<K)  =  SC*<CHNR<K>  -  MIR)  +  50 
CHMC(K)  =  SC:*(CHNC(K>  -  MIR)  +  50 
•  CONTINUE 

PERFORM  NEAREST  NEIGHBOR  INTERPOLATION 

DO  41  K  *  1»N 
ADDIR  =  0 
AJDOIC  »  0 

IF ( INT ( CHNR ( K ) ) ,  LT » INT ( CHNR ( K ) +0 ♦ 5 ) )  ADDIR 
IF  < INT ( CHNC  < K ) ) . LT , INT ( CHNC ( K ) +0 . 5 ) )  ADDIC 
CHNR < K >  "  INT (CHNR (K))  +  ADDIR 
CHNC(K)  ~  INT ( CHNC (K))  +  ADDIC 
CONTINUE 

WRITE  THE  CONTOUR  TO  IIS 


DO  93?  K  ■•=  1  .*  C  N-l  > 

1X1  =  INT < CHNR (K)) 

IY1  ~  INT ( CHNC < K ) ) 

1X2  =  INT ( CHNR ( K+l ) ) 

IY2  ■■=  INT  ( CHNC  (K+l  >  ) 

CALL  DVECT ( FCB ? 1X1 ? IY1 ? 1X2 1 IY2  ? NCHN ? NMASK ? BUFF ) 
CONTINUE 

.1X1  =  INT ( CHNR  CN)  > 

IY1  =  INT ( CHNC  <  M ) ) 

1X2  =  INT ( CHNR ( 1 ) ) 

IY2  =  INT ( CHNC ( 1 )  ) 

CALL  DVECT  <  FCB  ?  .1X1  ?  IY1  ?  .1X2  ?  I Y2  ?  NCHN  ?  NMASK BUFF ) 

TYPE  *? '  LEAVING  FD ' 

RETURN 

FORMAT (1 IX? 'DO  YOU  WANT  GRAPHICS  OR  REFRESH  DISPLAY 
FORMAT <1 IX? 'WHAT  COLOR  ?  ( 1 =RED ?  2=GREEN ?  4~BLUE)  '?*> 
FORMAT ( 1 IX ?  'DO  YOU  WANT  DISPLAY  ?  '  ?  >6 ) 

FORMAT <A2) 

END 


Atmospheric  Turbulence 

In  this  section,  we  will  examine  our  recent  work  in  the  area 
of  images  degraded  by  the  turbulent  atmosphere.  Our  efforts  have 
been  concentrated  in  the  areas  of  simulation  of  an  atmospheric 
turbulence  point  spread  function,  restoration  of  turbulence  de¬ 
grad  d  images  by  the  shi f t-and-add  method,  deblurring  of  shift-and- 
add  images,  and  the  effects  of  atmospheric  turbulence  on  the 
invariant  moments.  It  will  be  seen  that  the  combination  of  shift- 
and-add  processing  deblurring  substantially  improves  turbulence 
images. 

Background 

The  turbulent  atmosphere  degrades  images  of  objects  which  are 
observed  through  it.  These  degradations  are  due  to  the  passage 
of  the  optical  wavefront  through  regions  of  the  atmosphere  having 
variable  indices  of  refraction  caused  by  inhomogeneities  in 
temperature,  density,  and  velocity.  The  variable  indices  of  re¬ 
fraction  introduce  phase  and  amplitude  shifts  in  the  wavefront, 
although  for  relatively  mild  turbulence,  the  phase  shifts  are  the 
dominent  factor  in  image  degradation.  Phase  shifts  have  the 
effect  of  distributing  a  point  source  such  as  a  star  into  a 
speckle  image  as  in  Figure  II. 1. 

We  now  consider  an  object  f(x,y),  which  may  or  may  not  be  a 
point  source,  and  assume  a  mathematical  model  for  turbulence 
degradati  ons 
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(a)  Imaging  in  the  absence  of  turbulence  degradations. 


Point 

Source 

Object 


Turbulence 


Image 

(speckle  pattern) 


(b)  Imaging  tnrough  the  turbulent  atmosphere. 


Figure  1 


Sm(x,y)  *  hm(x,y)**f(x,y)  +  Cm(x,y) 

where  Sm(x,y)  is  the  mth  turbulence  image,  hm(x,y)  is  the  mth 
turbulence  point  spread  function  (speckle  pattern),  and  Cm(x,y) 
is  the  mth  contamination  term  encompassing  all  degradations 
such  as  system  nonlinearities  and  noise.  This  model  implicitly 
assumes  an  exposure  sufficiently  short  that  the  atmosphere  may 
be  considered  to  be  "frozen".  It  also  assumes  that  the  entire 
object  is  imaged  through  a  single  isoplanatic  patch  (region 
where  the  point  spread  function  may  be  considered  to  be  shift 
invariant).  The  next  section  presents  our  simulation  of  this 
model . 

Simulation 

In  order  to  carry  out  the  restoration  and  deblurring  of 
atmospheric  turbulence  images,  we  first  simulated  the  degraded 
images.  In  particular,  we  simulated  the  point  spread  functions 
hm  of  equation  (1)  which  were  then  convolved  with  our  original 
undegraded  Images.  At  this  point,  no  noise  was  added.  Our 
simulation  of  the  psf's  hm  followed  an  algorithm  developed  by 
B.  L.  McLamery  [1].  This  particular  algorithm  assumes  that  the 
phase  variations  produced  by  the  atmosphere  are  the  dominant 
components  of  image  degradation,  and  ignores  any  effects  due  to 
amplitude  variations.  We  summarize  this  method  as  follows: 

(1)  Generate  a  complex  array  of  Gaussian  random  numbers. 

(2)  Multiply  this  array  by  f"1^6,  the  square  root  of  the 


(1) 
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Kolmogorov  spatial  power  spectrum  (which  is  typically 
used  for  atmospheric  turbulence). 

(3)  Fourier  transform  the  array 

(4)  Separate  the  complex  result  of  this  transform  into 

its  real  and  imaginary  parts,  each  which  is  considered 
an  instaneous  phase  map  denoted  0(u,v). 

(5)  Choose  a  phase  map  0(u,v)  (from  step  4)  and  compute  the 


path  length  difference  map 


£(u,v)  = 


('  u ,  v )  A 
2  TT 


where  A  denotes  wavelength. 

(6)  Compute  the  coherent  point  spread  function 

s(x,y)  =  F  1  p(u, v)e 

where  we  have  at  this  point  set  d(u.v)s1,  and  F{ - >  denotes  Fourier  transform. 

(7)  Finally  compute  the  intensity  point  spread  function 

S(x,y)  as 

S(x,y)  =  |s(x,y)l2. 

This  S(x,y)  is  our  hm  of  equation  (1). 

Using  this  method,  we  started  at  step  (1)  with  arrays  of 
Gaussian  random  numbers  having  standard  deviations  of  0.5,  1.5, 

2.5,  and  3.5  and  generated  a  series  of  20  psf's  for  each  standard 
deviation.  Each  series  of  psf's  was  then  convolved  with  an  unde¬ 
graded  Image  to  produce  a  series  of  atmospheric  turbulence 
Images  on  which  the  shi f t-and-add  process  could  be  tested.  For 
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our  original  undegraded  images,  we  chose  a  drone  against  a  black 
(gray  level  0)  background.  We  chose  three  different  gray  levels  for 
the  drone:  64,  128  and  192.  Also,  because  the  shi f t-and-add 
processing  requires  a  "brightest  point"  in  the  original  object, 
a  "glint"  or  point  source  at  a  gray  level  of  255  was  inserted  on 
the  drone  body. 

Results  of  this  simulation  are  depicted  in  Figures  2-4. 

Figure  2  shows  the  original  drone  at  a  gray  level  of  128.  Figure 
3  shows  examples  from  the  series  of  20  atmospheric  turbulence 
psf's.  Finally,  Figure  4  depicts  the  corresponding  examples 
from  the  series  of  turbulence  degraded  images. 

Restoration  by  Shi ft-and-Add  Processing 

Having  simulated  atmospheric  turbulence  images,  we  attempted 
restoration  by  the  sh i f t-and-add  method  developed  by  Bates  and 
Cady  [2,3].  We  continue  to  assume  the  degradation  model  of 
equation  (1) 

S  (x,y)  3  h  (x,y)**f(x,y)  +  C  (x,y) 
mm  m 

where,  as  before,  S  (x,y)  denotes  the  mth  turbulence  image.  We 

m 

also  assume  that  we  have  M  of  these  images.  The  shi  fv  and-add 

process  is  then  simple  to  describe.  Let  (£  ,y  )  denote  the 

m  m 

spatial  coordinates  at  which  the  degraded  image  S  (x,y)  takes  on 

its  maximum  value.  The  shi f t-and-add  image  Is  the  result  of 

translating  each  S  (x,y)  so  that  U  »y  )  is  at  the  origin  of 

m  m  m 

coordinates  and  summing  the  M  translated  Images: 

*  i  ?  S_( x  +  £_,y  +  U  ). 


s(x,y) 


f  o\ 


3(a)  Original  s.d.  of  Gaussian  random  numbers  was  1.5 


- 


3(b)  Original  s.d.  of  Gaussian  random  numbers  was  2.5 


4(a)  Original  s.d.  of  Gaussian  random  numbers  was  1.5 


4(b)  Original  s.d.  of  Gaussian  random  numbers  was  2.5 


4(c)  Original  s.d.  of  Gaussian  random  numbers  was  3.5. 


Figure  4:  Simulated  turbulence  images  (magnified  4 x ) . 


Because  all  operations  involved  are  linear,  we  may  also  write  the 
shi f t-and-add  image  s{x,y)  as 


s(x,y)  =  h(x,y)**f(x,y)  +  c(x,y)  (3) 

where  h(x,y)  is  the  sum  of  the  shifted  hm(x,y)  and  c(x,y)  is 
the  sum  of  the  shifted  c(x,y). 

A  mathematical  analysis  of  this  method  has  been  performed 
by  Hunt,  Fright,  and  Bates  [4].  This  analysis  indicates  that  for 
a  large  number  M  of  speckle  images,  a  spike  corresponding  to  the 
"brightest  point"  in  the  original  object  will  be  created  near  the 
origin  of  coordinates.  The  proper  alignment  of  the  other  portions 
of  the  image  is  due  to  the  linear  nature  of  convolution.  This 
process  may  break  down  when  a  "brightest  point"  is  unavailable 
in  the  original  object  or  when  several  points  or  areas  of  ap¬ 
proximately  equal  brightness  are  present.  In  this  case,  the  de¬ 
graded  images  may  be  incorrectly  aligned  causing  "ghosts"  to 
appear  in  the  result. 

We  performed  shi f t-and-add  processing  on  each  series  of  20 
degraded  images  which  had  been  simulated  by  the  method  of  section 
II. 2,  and  the  results  are  summarized  in  Table  1.  We  see  from 
this  table  that  both  the  contrast  between  the  point  source  (or 
"glint")  and  the  aircraft  body  and  also  the  original  standard 
deviation  of  the  Gaussian  random  numbers  used  in  computing  the 
turbulence  psf  affect  the  performance  of  shi f t-and-add  processing. 
It  does  appear  however,  that  even  in  the  cases  where  shift-and- 
add  processing  produced  a  blurred  resultant  image,  the  shift-and- 
add  image  should  be  a  better  candidate  for  deblurring  than  are 
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reconstruction 
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reconstruct! on 

3.5 

highly  blurred 
image  reconstruc¬ 
tion 

highly  blurred 
image  reconstruc- 
ti  on 

highly  blurred 
image  reconstruc- 
ti  on 

Table  1 


the  individual  turbulence  images.  Shi f t-and-add  images  for  the 
case  of  an  original  drone  gray  level  of  128  are  depicted  in  Figure  5. 

We  also  ran  one  test  to  check  the  effect  of  increasing  the 
number  M  of  processed  images  M  from  20  to  50.  The  case  chosen 
was  that  of  the  aircraft  gray  level  at  128  and  the  turbulence  psf 
generated  from  Gaussian  random  numbers  with  a  standard  deviation 
of  3.5.  It  is  interesting  to  note  that  for  this  case,  the  shift- 
and-add  images  for  M  =  20  and  M  3  50  are  visually  indistinguishable. 

It  would  therefore  appear  likely  that  beyond  some  point,  for  this 
case  M  =  20  or  perhaps  less,  no  improvement  is  gained  by  processing 
more  images. 

Debl urri ng 

As  noted  in  the  previous  section,  the  results  of  shift-and- 
add  processing  tended  to  be  blurred,  especially  for  higher  levels 
of  turbulence.  Assuming  the  model  of  equation  (3)  for  the  shift- 
and-add  image 

s(x,y)  =  h(x,y)**f(x,y)  +  c(x,y), 
we  wish  to  recover  f(x,y).  Our  task  is  somewhat  simplified  by 
the  fact  that  our  simulated  Images  are  free  of  degradations  other 
than  atmospheric  turbulence,  i.e.,  we  may  essentially  ignore  the 
term  c(x,y).  Thus  the  removal  of  blur  involves  estimation  of 
the  cumulative  point  spread  function  h(x,y)  and  filtering. 

To  estimate  the  cumulative  psf,  we  inserted  a  point  source  < 

near  but  outside  the  drone  body  (Figure  6).  This  point  source 
was  at  a  lower  gray  level  than  the  one  of  the  drone  body  so  that 
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it  would  not  interfere  with  the  shi f t-and-add  processing.  Series 
of  20  turbulence  images  were  then  created  by  convolution  with  the 
simulated  point  spread  functions  hm(x,y),  and  shi ft-and-add  images 
were  formed  from  each  of  these  series.  By  checking  the  results  of 
this  processing  on  the  second  point  source,  we  were  able  to  obtain 
a  point  spread  function  h(x,y)  for  the  entire  process.  Using  a 
least  squares  fit,  this  point  spread  function  was  approximated  by 
a  Gaussian  to  be  used  in  filtering. 

In  order  to  recover  the  original  object,  we  used  the  Wiener 
and  Cannon  filters.  The  form  of  the  Wiener  filter  used  was 


W  ( u ,  v ) 


H* ( u , v ) 

|H(u,v)|2  *  an2 

*fiu7v) 


where  H(u,v)  is  the  Fourier  transform  of  the  Gaussian  approxima¬ 
tion  to  h(x,y),  *  denotes  complex  conjugation,  |'|  denotes  the 

p 

modulus  of  the  transform,  on  is  noise  variance,  and  4>f(u,v)  is 
the  ideal  image  power  spectrum.  The  results  of  this  filtering 
are  shown  in  Figure  7,  and  we  see  a  great  improvement  in  sharpness 
of  the  drone  outline,  especially  for  the  cases  of  2.5  and  3.5 
standard  deviation  in  the  original  Gaussian  iterates.  We  also 
note  fairly  severe  "ringing"  in  these  filtered  images  -  more  than 
Is  actually  visible  in  Figure  7.  Because  we  are  particualrly 
Interested  in  obtaining  the  silhouette  or  outline  of  the  object 
for  use  in  calculating  moments  or  Fourier  descriptors.  It  is 
advantageous  to  accept  some  ringing  In  order  to  retain  sharpness 
of  the  outline.  We  thresholded  the  restored  Images  in  order  to 
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7(c)  Original  s.d.  of  Gaussian  random  numbers  is  3.5. 


Figure  7 


Wiener  filter  deblurring  of  shi f t-and-add  images 
(magni fied  4x) . 


compare  silhouettes  (or  outlines)  with  that  of  the  thresholded 
original.  These  images  are  shown  in  Figures  8  and  9. 

We  also  tried  deblurring  with  the  Cannon  filter  which  some 
times  increases  object  sharpness.  We  used  the  form 

C( u,v) 

In  our  case,  the  results  from  this  filter  were  poorer  than  those 
from  the  Wiener  filter. 

From  Figures  4  and  9,  we  can  see  that  the  combination  of 
shi f t-and-add  processing  and  Wiener  filter  deblurring  produces 
marked  improvements  in  images  degraded  by  atmospheric  turbulence. 

Moments 

We  have  previously  run  te<ts  to  determine  the  effect  of  rota¬ 
tion,  scale  change,  and  noise  on  invariant  moments.  These  results 
indicated  that  any  classification  system  based  on  moments  should 
weight  the  third  and  seventh  moments  less  heavily  and  should  use 
the  largest  weights  on  the  first  and  second  moments. 

We  have  now  considered  the  effect  of  atmospheric  turbulence 
on  the  invariant  moments.  Our  first  test  involved  the  effect 
of  atmospheric  turbulence  and  of  shi f t-and-add  processing  alone. 
For  each  standard  deviation,  moments  were  calculated  for  the 
original  undegraded  image  (aircraft  gray  level  128),  for  two  of 
the  turbulence  Images  from  the  series  of  20,  and  for  the  re¬ 
sulting  shi ft-and-add  image.  The  percent  error  was  measured 


Figure  8:  Thresholded  original  drone  (magnified  4 x ) . 


9(a)  Original  s.d.  of  Gaussian  random  numbers  was  1.5 


9(b)  Original  s.d.  of  Gaussian  random  numbers  was  2.5 


9(c)  Original  s.d.  of  Gaussian  random  numbers  was  3.5. 


Figure  9:  Thresholded  restored  images  (magnified  4x). 


relative  to  the  moments  of  the  original  image  and  is  plotted  in 
Figures  10-13.  From  these  figures,  it  appears  that  the  first  and 
especially  the  fifth  moment  are  most  sensitive  to  turbulence  de¬ 
gradations  and  possibly  should  be  weighted  less  heavily  in  any 
classification  schemes.  These  figures  also  show  that  contrary 
to  expectation,  an  original  standard  deviation  of  1.5  produces 
the  most  error  in  the  moments,  more  than  the  more  severe  degrada¬ 
tions  produced  by  standard  deviations  of  2.5  and  3.5.  The  moments 
of  the  shi f t-and-add  image  for  this  case  are  also  higher  in  error 
than  are  those  of  the  other  shi f t-and-add  images.  We  may  explain 
this  by  referring  to  Table  1  where  we  note  the  presence  of  a  faint 
"ghost"  aircraft  in  the  background.  We  do  note  that  in  all  cases > 
except  for  1.5  standard  deviation  in  the  original  Gaussian  iterates 
shi f t-and-add  processing  did  reduce  the  error  in  the  moments. 

Our  second  test  involved  calculating  the  moments  of  the  binary 
thresholded  original  image  (Figure  8)  and  of  the  shi  f t-and-add 
images  after  Wiener  filtering  and  thresholding  (Figure  9).  This 
added  processing  actually  brought  about  little  reduction  in  the 
error  of  the  moments  as  may  be  seen  in  Figure  14.  However,  the 
accuracy  of  the  first  moment  was  i  mproved,  whi  ch  is  fortunate, 
since  our  previous  tests  involving  rotation,  size  change,  and  noise 
indicated  that  the  first  moment  should  be  weighted  heavily  in 
classification  algorithms.  As  in  the  case  for  turbulence  or  shift- 
and-add  images,  the  fifth  moment  appears  to  be  most  subject  to 
error  and  should  thus  be  weighted  less  heavily  in  any  classifica¬ 
tion  schemes  based  on  moments. 
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| %  ERROR |  vs.  MOMENT 
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Original  s.d.  of  the  Gaussian  Iterates  used  In  computing  the 
turbulence  psf's  Is  3.5 

-  turbulence  corrupted  Image  #5 

•  turbulence  corrupted  Image  #15 

shl ft-and-add  Image  (20  files) 
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relative  to  the  moments 
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